version 11
#delim ;
set more off;
capture log close;
capture clear;
local adodir "C:\Scott\Synch\Stata\ado";
sysdir set PERSONAL "`adodir'\personal";
sysdir set PLUS "`adodir'\plus";
sysdir set OLDPLACE "`adodir'";


/***!***!***!***!***!*** [National_bank_analysis2.do ] ***!***!***!***!***!
*
* Project: National Banks 		
* Programmer:  Scott Fulford
*
* Date:    	 7/12/2010
*
* Auditor:      
* Audit Date:   
*
* Purpose:      
* 1) Create a multiple imputation data set based 1890 county level data of census and National Banking which imputes the optimal level of banking if the capital was unconstrained
* 2) 
* 3)
* Inputs: 
      National_Banks_counties1890_addvar created by National_Banks_analysis1.do
*
* Ouputs: 
*		MI_National_banks
*
*
***!***!***!***!***!***!***!***!***!***!***!***!***!***!***!***!***/


/***Define Global Directory ****/
	local INDIR "C:\Scott\Research\National_Banks\Intermediate";
	local PROGDIR  "C:\Scott\Research\National_Banks\Programs";
	local OUTDIR  "C:\Scott\Research\National_Banks\Intermediate";
	local GRAPHDIR "C:\Scott\Research\National_Banks\Intermediate";
	local TEXDIR "C:\Scott\Research\National_Banks\Intermediate\Regression_output";

/***Define Sample selection***/
/*
Select the sample by saving different MI files.
allrural selects all rural counties (including southern and western ones)
unionrural selects only "union" rural counties (the base regressions)
*/
/*local sampleselect allrural;	*/
local sampleselect unionrural;	
/*******************************/

*log using National_bank_analysis2_`sampleselect'.txt, text replace;

set memory 800m;
set seed 106916162; /*Seed generated from 1-500,000,000 by www.random.org from atmospheric noise*/

/***** Create replications file *****/
local replications = 100;
local Ninitial = 5000; /*How many steps between draws in Gibbs Sampler*/

/***** When to end search for best capital stock */
local tolerance 10^-3;

cd `INDIR';
use National_Banks_counties1890_addvar, clear;
capture drop __000001; /*Stata has some sort of problem not dropping a tempvar*/

/*********
Create estimate of under/over capitalization
**********/
keep if insample_`sampleselect';
/*
set trace on;
set tracedepth 1;
*/

cap program drop transitioncap;
program transitioncap, rclass;
	version 11;
	args todo idivide ;
	tempname lnf;

	if `idivide' >=50 | `idivide' <= 0  {;
		scalar `lnf' = .;
	}; 
	else {;

	
		tempvar lncp_up lncp_down;
		/*Upper and lower bounds if exact capitalization*/
		quietly {;
		gen `lncp_up' = ln_capitalstock_pc;
		gen `lncp_down' = `lncp_up'; /*Avoid any last digit errors*/	
		
		/*Lower bound if have cap 50 banks*/
		replace `lncp_down' = ln( (capitalstock - (50-`idivide')*banks50c)/tpop) if banks50c>0 & banks50c <.; 
		/*Upper and lower bound if missing*/
		replace `lncp_up' = ln(`idivide'/tpop) if banks==0;
		replace `lncp_down'=. if banks==0;
		};
		xtintreg `lncp_down' `lncp_up' i.year ln_tpop, intreg;
		/*Report log likelihood*/	
		scalar `lnf' = `e(ll)'; /*Only report a single value since it is the overall log likelihood and no sum is necessary*/
	};
	return scalar lnf = `lnf';
end;

cap program drop transitioncap25;
program transitioncap25, rclass;
	version 11;
	args todo idivide idivide25 ;
	tempname lnf;

	if (`idivide' >=50 | `idivide' <= 0) | (`idivide25' >=25 | `idivide' <= 0)   {;
		scalar `lnf' = .;
	}; 
	else {;
		tempvar lncp_up lncp_down;
		quietly {;
			
		/*Upper and lower bounds if exact capitalization*/
			
		gen `lncp_up' = ln_capitalstock_pc;
		replace `lncp_up'= ln(capitalstock1902/tpop) if year ==1900 ;
		gen `lncp_down' = `lncp_up'; /*Avoid any last digit errors*/
		
		/*Lower bound if have cap 50 banks*/
		replace `lncp_down' = ln( (capitalstock - (50-`idivide')*banks50c)/tpop) if banks50c>0 & banks50c <. & year !=1900; 
		/*Lower bound if have cap 25 banks in 1900*/
		replace `lncp_down' = ln( (capitalstock1902 - (25-`idivide25')*banks25c1902)/tpop) if banks25c1902>0 & banks25c1902 <. & year ==1900; 
		
		/*Upper and lower bound if missing*/
		replace `lncp_up' = ln(`idivide'/tpop) if banks==0 & year !=1900;
		replace `lncp_up' = ln(`idivide25'/tpop) if banks1902==0 & year ==1900 ;
		
		replace `lncp_down'=. if banks==0 & year !=1900 ;
		replace `lncp_down'=. if banks1902==0 & year ==1900 ;
		};
		xtintreg `lncp_down' `lncp_up' i.year ln_tpop, intreg;
		/*Report log likelihood*/	
		scalar `lnf' = `e(ll)'; /*Only report a single value since it is the overall log likelihood and no sum is necessary*/
	};
	return scalar lnf = `lnf';
end;


/*************
Estimate Maximum Likelihood 
	and find best transistions
*************/

/***1902*****/
cd `INDIR';
use National_Banks_counties1890_addvar, clear;
capture drop __000001; /*Stata has some sort of problem not dropping a tempvar*/
/*********
Create estimate of under/over capitalization with 25 banks
**********/

keep if insample_`sampleselect' & capitalstock<=500;

local bestloglikelihood = -10^100; /*A big number*/

local adivide1 = 12;
local bdivide1 = 30;
local adivide2 = 5;
local bdivide2 = 20;


local steps = 10;
local stepsize1 =(`bdivide1'-`adivide1')/`steps';
local stepsize2 =(`bdivide2'-`adivide2')/`steps';


while (max(`stepsize1',`stepsize2')>`tolerance') {;
	local stepsize1 =(`bdivide1'-`adivide1')/`steps';
	local stepsize2 =(`bdivide2'-`adivide2')/`steps';
	forvalues idivide1 = `adivide1'(`stepsize1')`bdivide1' {;
		local bdivide2_hat = min(`idivide1',`bdivide2'); /*Dividing line for 25 capital should be less than 50 capital, enforce that and reduce search*/
		
		forvalues idivide2 = `adivide2'(`stepsize2')`bdivide2_hat' {;
			
			quietly transitioncap25 0 `idivide1' `idivide2';
			if `r(lnf)' >= `bestloglikelihood' {;
				local bestloglikelihood = `r(lnf)';
				local bestdivide1 = `idivide1';
				local bestdivide2 = `idivide2';
			};
			noisily disp "Divide50=" `idivide1' ", Divide25=" `idivide2' ",  Likelihood=" `r(lnf)';
		};
	};
	local adivide1 = `bestdivide1'- `stepsize1';
	local bdivide1 = `bestdivide1' +`stepsize1';
	
	local adivide2 = `bestdivide2'- `stepsize2';
	local bdivide2 = `bestdivide2' +`stepsize2';
};
disp "Best divide50=" `bestdivide1';
disp "Best divide25=" `bestdivide2';
disp "Best LL=" `bestloglikelihood';

/***Do final again to get coefficients*/
transitioncap25 0 `bestdivide1' `bestdivide2';
local sigma25_eta= e(sigma_u);
local sigma25_eps= e(sigma_e);
matrix eb= e(b);
local alpha25 = eb[1,colnumb(eb,"_cons")];
local alpha25_1880 = eb[1,colnumb(eb,"1880.year")];
local alpha25_1890 = eb[1,colnumb(eb,"1890.year")];
local alpha25_1900 = eb[1,colnumb(eb,"1900.year")];
local bestdivide25 = `bestdivide2';
local bestdivide50 = `bestdivide1';
local theta25 = eb[1,colnumb(eb,"ln_tpop")];

/*quadchk;*/
cd "`TEXDIR'";
outreg2 using excess_capital25_`sampleselect', excel replace;


/*********
Create estimate of under/over capitalization for 1870-1900 without 25 banks
**********/

local bestloglikelihood = -10^100; /*A big number*/
local initialdivide =20;
local adivide = 12;
local bdivide = 30;
local steps = 10;
local stepsize =(`bdivide'-`adivide')/`steps';

while (`stepsize'>`tolerance') {;
	local stepsize =(`bdivide'-`adivide')/`steps';	
	forvalues idivide = `adivide'(`stepsize')`bdivide' {;
		quietly transitioncap 0 `idivide';
		if `r(lnf)' >= `bestloglikelihood' {;
			local bestloglikelihood = `r(lnf)';
			local bestdivide = `idivide';
		};
		noisily disp "Divide=" `idivide' ",  Likelihood=" `r(lnf)';
	};
	local adivide = `bestdivide'- `stepsize';
	local bdivide = `bestdivide' +`stepsize';
};
disp "Best divide=" `bestdivide';
disp "Best LL=" `bestloglikelihood';

/***Do final again to get coefficients*/
transitioncap 0 `bestdivide';
local sigma_eta= e(sigma_u);
local sigma_eps= e(sigma_e);
matrix eb= e(b);
local alpha = eb[1,colnumb(eb,"_cons")];
local alpha_1880 = eb[1,colnumb(eb,"1880.year")];
local alpha_1890 = eb[1,colnumb(eb,"1890.year")];
local alpha_1900 = eb[1,colnumb(eb,"1900.year")];
local theta = eb[1,colnumb(eb,"ln_tpop")];
/*quadchk;*/
cd "`TEXDIR'";
outreg2 using excess_capital_`sampleselect', excel replace;


/*********
Old code to find optimal parameter
*****

/*Do bisection search to find optimal parameter*/
local tolerance = 10^-3;

local adivide = `bestdivide'- 1;
local bdivide = `bestdivide' + 1;
local alpha1 = (3-sqrt(5))/2;
local alpha2 = (sqrt(5)-1)/2;
local x1 = `adivide' + `alpha1'*(`bdivide'-`adivide');
	transitioncap 0 `x1';
	local f1 =`r(lnf)';
local x2 = `adivide' + `alpha2'*(`bdivide'-`adivide');
	transitioncap 0 `x2';
	local f2 =`r(lnf)';
local diff = `alpha1'*`alpha2'*(`bdivide'-`adivide');
noisily disp "Divide1= `x1' Divide1= `x2'  ";
while   (`diff' >`tolerance') {;
	local diff = `diff'*`alpha2';
	if (`f2'<`f1') {;
			local x2 = `x1';
			local x1 = `x1' - `diff';
			local f2 = `f1';
			transitioncap 0 `x1';
			local f1 =`r(lnf)';
		}; else {;
			local x1 = `x2';
			local x2 = `x2' +`diff';
			local f1 = `f2';
			transitioncap 0 `x2';
			local f2 =`r(lnf)';
		};
		noisily disp "Divide1= `x1' Divide1= `x2'  ";
};
if (`f2'<`f1') {;
		local bestdivide = `x2';
	}; else {;
		local bestdivide = `x1';
	};
*****/
/* ***For testing *****
/****For drawing capital ******/
if "`sampleselect'" == "allrural" {;
/*only 50*/
local bestdivide = 17.041;
local sigma_eta= 1.76;
local sigma_eps= 0.7197;
local alpha = -15.239;
local alpha_1880 = -.0024;
local alpha_1890 = 0.5742;
local alpha_1900 = 0.2702;
local theta = .8098;
/*25 in 1902*/
local bestdivide25 = 13.2;
local bestdivide50 =22.77;
local sigma25_eta= 1.68;
local sigma25_eps= 0.709;
local alpha25 = -12.94;
local alpha25_1880 = 0.0045;
local alpha25_1890 = 0.5327;
local alpha25_1900 = 0.5349;
local theta25 = .6018;

}; else if "`sampleselect'" == "unionrural" {;
/*only 50*/
local bestdivide = 17.041;
local sigma_eta= 1.76;
local sigma_eps= 0.7197;
local alpha = -15.239;
local alpha_1880 = -.0024;
local alpha_1890 = 0.5742;
local alpha_1900 = 0.2702;
local theta = .8098;
/*25 in 1902*/
local bestdivide25 = 13.2;
local bestdivide50 =22.77;
local sigma25_eta= 1.68;
local sigma25_eps= 0.709;
local alpha25 = -12.94;
local alpha25_1880 = 0.0045;
local alpha25_1890 = 0.5327;
local alpha25_1900 = 0.5349;
local theta25 = .6018;

}; else {;
	error "Not recognized sample";
};
*/


/******** 
Using log likelihood parameters, create bounds for capital
Which counties have too little, too much capital?

********/

cd `INDIR';
use National_Banks_counties1890_addvar, clear;
capture drop __000001; /*Stata has some sort of problem not dropping a tempvar*/
keep if insample_`sampleselect';
xi i.year; /*Create year dummies*/

gen needscapital = ( (year >=1870 & year <=1900) &
				((banks50c>0 & banks50c<.)/*Counties with a 50 bank*/
				| (banks==0) ))  /*Counties without a bank*/
				if banks <.; 



gen lowcapital = 
		/*Counties with a 50 bank*/
			ln(capitalstock - (50-`bestdivide')*banks50c)
			-ln(tpop)-`alpha' 
						- _Iyear_1880*`alpha_1880'
						- _Iyear_1890*`alpha_1890'
						- _Iyear_1900*`alpha_1900'
						- `theta'*ln_tpop
				if (banks50c>0 & banks50c<.);
	gen highcapital =
			/*Counties with a 50 bank*/
					ln(capitalstock)
					-ln(tpop)-`alpha' 
							- _Iyear_1880*`alpha_1880'
							- _Iyear_1890*`alpha_1890'
							- _Iyear_1900*`alpha_1900'
							- `theta'*ln_tpop
					if (banks50c>0 & banks50c<.);
replace highcapital = 
			/*Counties with a no banks*/
				ln(`bestdivide')
				-ln(tpop)-`alpha' 
						- _Iyear_1880*`alpha_1880'
						- _Iyear_1890*`alpha_1890'
						- _Iyear_1900*`alpha_1900'
						- `theta'*ln_tpop
				if (banks==0); 				

/** For 1902 data with 25. Note that only 1900 has 25 capital banks, but have a new
estimate of the divide for 1870, 1880 and 1890, so a new set of dividers*/
gen needscapital25 = needscapital; /*If only change for those counties which do 
not have a bank in 1900 but get one in 1902, then same counties need a capital draw, just from a different distribution */
/*replace needscapital25 = (
				((banks25c1902>0 & banks25c1902<.)/*Counties with a 25 bank*/
				| (banks1902==0) ))  /*Counties without a bank*/
				if banks <. & year == 1900; */

/* For 1870, 1880, 1890, and 1900 for counties without 25 bank*/
gen lowcapital25 = 
		/*Counties with a 50 bank*/
				ln((capitalstock - (50-`bestdivide50')*banks50c))
				-ln(tpop)-`alpha25' 
						- _Iyear_1880*`alpha25_1880'
						- _Iyear_1890*`alpha25_1890'
						- _Iyear_1900*`alpha25_1900'
						- `theta25'*ln_tpop
				if (banks50c>0 & banks50c<.);
gen highcapital25 =
		/*Counties with a 50 bank*/
				ln(capitalstock)  
				-ln(tpop)-`alpha25' 
						- _Iyear_1880*`alpha25_1880'
						- _Iyear_1890*`alpha25_1890'
						- _Iyear_1900*`alpha25_1900'
						- `theta25'*ln_tpop
				if (banks50c>0 & banks50c<.);
replace highcapital25 = 
			/*Counties with a no banks*/
				ln(`bestdivide50')
				-ln(tpop)-`alpha25' 
						- _Iyear_1880*`alpha25_1880'
						- _Iyear_1890*`alpha25_1890'
						- _Iyear_1900*`alpha25_1900'
						- `theta25'*ln_tpop
				if (banks==0); 			
				
/*Change 1902 values for 1900, and adjust for 25 capital banks
Note that only change for those that had only one 25 capital bank*/
replace lowcapital25 = 
		/*Counties with a 25 bank*/
				ln(`bestdivide25')
				-ln(tpop)-`alpha25' 
						- _Iyear_1880*`alpha25_1880'
						- _Iyear_1890*`alpha25_1890'
						- _Iyear_1900*`alpha25_1900'
						- `theta25'*ln_tpop
						 if (banks25c1902==1 & banks==0 /*Banks is 1900 banks*/) & year == 1900;
replace highcapital25 =
		/*Counties with a 25 bank in 1902 could have entered with a 50 bank in 1900
		so highest capital is min of 25 and bestdivide50*/
				ln(min(capitalstock1902,`bestdivide50' )) 
				-ln(tpop)-`alpha25' 
						- _Iyear_1880*`alpha25_1880'
						- _Iyear_1890*`alpha25_1890'
						- _Iyear_1900*`alpha25_1900'
						- `theta25'*ln_tpop
				if (banks25c1902==1 & banks==0)  & year == 1900;
replace highcapital25 = 
			/*Counties with a no banks*/
				ln(`bestdivide25')
				-ln(tpop)-`alpha25' 
						- _Iyear_1880*`alpha25_1880'
						- _Iyear_1890*`alpha25_1890'
						- _Iyear_1900*`alpha25_1900'
						- `theta25'*ln_tpop
				if (banks1902==0 & banks==0) & year == 1900; 	 			


/*

/*by gisjoin1890: egen stillneedscapital = max(needscapital); Will be missing if any missing in county, that should be correct behavior*/


gen lowcapital = 
		/*Counties with a 50 bank*/
				ln((capitalstock - (50-`bestdivide')*banks50c)/tpop)-`alpha' 
						- _Iyear_1880*`alpha_1880'
						- _Iyear_1890*`alpha_1890'
						- _Iyear_1900*`alpha_1900'
				if (banks50c>0 & banks50c<.);
gen highcapital =
		/*Counties with a 50 bank*/
				ln(capitalstock /tpop)-`alpha' 
						- _Iyear_1880*`alpha_1880'
						- _Iyear_1890*`alpha_1890'
						- _Iyear_1900*`alpha_1900'
				if (banks50c>0 & banks50c<.);
replace highcapital = 
			/*Counties with a no banks*/
				ln(`bestdivide'/tpop)-`alpha' 
						- _Iyear_1880*`alpha_1880'
						- _Iyear_1890*`alpha_1890'
						- _Iyear_1900*`alpha_1900'
				if (banks==0); 				

/** For 1902 data with 25. Note that only 1900 has 25 capital banks, but have a new
estimate of the divide for 1870, 1880 and 1890, so a new set of dividers*/
gen needscapital25 = needscapital;
replace needscapital25 = (
				((banks25c1902>0 & banks25c1902<.)/*Counties with a 25 bank*/
				| (banks1902==0) ))  /*Counties without a bank*/
				if banks <. & year == 1900; 

/* For 1870, 1880, 1890*/
gen lowcapital25 = 
		/*Counties with a 50 bank*/
				ln((capitalstock - (50-`bestdivide50')*banks50c)/tpop)-`alpha25' 
						- _Iyear_1880*`alpha25_1880'
						- _Iyear_1890*`alpha25_1890'
						- _Iyear_1900*`alpha25_1900'
				if (banks50c>0 & banks50c<.);
gen highcapital25 =
		/*Counties with a 50 bank*/
				ln(capitalstock /tpop)-`alpha25' 
						- _Iyear_1880*`alpha25_1880'
						- _Iyear_1890*`alpha25_1890'
						- _Iyear_1900*`alpha25_1900'
				if (banks50c>0 & banks50c<.);
replace highcapital25 = 
			/*Counties with a no banks*/
				ln(`bestdivide50'/tpop)-`alpha25' 
						- _Iyear_1880*`alpha25_1880'
						- _Iyear_1890*`alpha25_1890'
						- _Iyear_1900*`alpha25_1900'
				if (banks==0); 			
				
/*Change 1902 values for 1900, and adjust for 25 capital banks*/
replace lowcapital25 = 
		/*Counties with a 25 bank*/
				ln((capitalstock1902 - (25-`bestdivide25')*banks25c1902)/tpop)  -`alpha25' 
						- _Iyear_1880*`alpha25_1880'
						- _Iyear_1890*`alpha25_1890'
						- _Iyear_1900*`alpha25_1900'
				if (banks25c1902>0 & banks25c1902<.) & year == 1900;
replace highcapital25 =
		/*Counties with a 25 bank*/
				ln(capitalstock1902  /tpop)-`alpha25' 
						- _Iyear_1880*`alpha25_1880'
						- _Iyear_1890*`alpha25_1890'
						- _Iyear_1900*`alpha25_1900'
				if (banks25c1902>0 & banks25c1902<.)  & year == 1900;
replace highcapital25 = 
			/*Counties with a no banks*/
				ln(`bestdivide25'/tpop)-`alpha25' 
						- _Iyear_1880*`alpha25_1880'
						- _Iyear_1890*`alpha25_1890'
						- _Iyear_1900*`alpha25_1900'
				if (banks1902==0) & year == 1900; 	 			
*/

/******** Generate unobserved values *********/

/*Stata and Mata don't get along perfectly, can't invoke mata in a loop, only issue single command*/
mata:;
	void makedraws(real scalar Ninitial) {
	/*Ninitial Number of interations to get to run before take draw*/
	
	
	replications = strtoreal(st_local("replications"));
	sigma_eps=strtoreal(st_local("sigma_eps"));
	sigma_eta=strtoreal(st_local("sigma_eta"));
	etai_c = st_data(., "etai_c");
  lowcapital1870=st_data(. ,"lowcapital1870");
  lowcapital1880=st_data(. ,"lowcapital1880");
  lowcapital1890=st_data(. ,"lowcapital1890");
  lowcapital1900=st_data(. ,"lowcapital1900");    
  highcapital1870=st_data(. ,"highcapital1870");
  highcapital1880=st_data(. ,"highcapital1880");
  highcapital1890=st_data(. ,"highcapital1890");
  highcapital1900=st_data(. ,"highcapital1900");
	
	numberdraws = rows(etai_c);
	/*Draw for each year*/
	epsi_c1870 =(rtnorm(1,1, J(1,numberdraws,0), J(1,numberdraws,sigma_eps),(lowcapital1870-etai_c)',(highcapital1870-etai_c)'))';
	epsi_c1880 =(rtnorm(1,1, J(1,numberdraws,0), J(1,numberdraws,sigma_eps),(lowcapital1880-etai_c)',(highcapital1880-etai_c)'))';
	epsi_c1890 =(rtnorm(1,1, J(1,numberdraws,0), J(1,numberdraws,sigma_eps),(lowcapital1890-etai_c)',(highcapital1890-etai_c)'))';
	epsi_c1900 =(rtnorm(1,1, J(1,numberdraws,0), J(1,numberdraws,sigma_eps),(lowcapital1900-etai_c)',(highcapital1900-etai_c)'))';
	
	for (repnum =1; repnum<=replications; repnum++) {
		for (xx=1; xx<=Ninitial; xx++) {
			lbeta=rowmax( (lowcapital1870-epsi_c1870, lowcapital1880-epsi_c1880, lowcapital1890-epsi_c1890, lowcapital1900-epsi_c1900));
		
			ubeta=rowmin( (highcapital1870-epsi_c1870, highcapital1880-epsi_c1880, highcapital1890-epsi_c1890, highcapital1900-epsi_c1900));
			
			etai_c =(rtnorm(1,1, J(1,numberdraws,0), J(1,numberdraws,sigma_eta), lbeta', ubeta'))';
		
			/*Draw for each year*/
			epsi_c1870 =(rtnorm(1,1, J(1,numberdraws,0), J(1,numberdraws,sigma_eps),(lowcapital1870-etai_c)',(highcapital1870-etai_c)'))';
			epsi_c1880 =(rtnorm(1,1, J(1,numberdraws,0), J(1,numberdraws,sigma_eps),(lowcapital1880-etai_c)',(highcapital1880-etai_c)'))';
			epsi_c1890 =(rtnorm(1,1, J(1,numberdraws,0), J(1,numberdraws,sigma_eps),(lowcapital1890-etai_c)',(highcapital1890-etai_c)'))';
			epsi_c1900 =(rtnorm(1,1, J(1,numberdraws,0), J(1,numberdraws,sigma_eps),(lowcapital1900-etai_c)',(highcapital1900-etai_c)'))';
		}
		st_store(.,("eps"+strofreal(repnum)+"_c1870", "eps"+strofreal(repnum)+"_c1880", "eps"+strofreal(repnum)+"_c1890", "eps"+strofreal(repnum)+"_c1900", "eta"+strofreal(repnum)+"_c"),(epsi_c1870, epsi_c1880, epsi_c1890, epsi_c1900, etai_c ) );
	};
}
end

/*Get starting value for eta, better if reasonable*/

gen eta0ct =lowcapital+ (highcapital-lowcapital)/2 if banks50c>0;
replace eta0ct = ln(`bestdivide'/(2*tpop))-`alpha' 
						- _Iyear_1880*`alpha_1880'
						- _Iyear_1890*`alpha_1890'
						- _Iyear_1900*`alpha_1900'
						- `theta'*ln_tpop
				if (banks==0);
by gisjoin1890: egen etai_c = mean(eta0ct);
drop eta0ct; 

gen eta0ct25 =lowcapital25+ (highcapital25-lowcapital25)/2 if banks50c>0;
replace eta0ct25 =lowcapital25+ (highcapital25-lowcapital25)/2 if banks25c1902>0 & year == 1900;
replace eta0ct25 = ln(`bestdivide50'/(2*tpop))-`alpha25' 
						- _Iyear_1880*`alpha25_1880'
						- _Iyear_1890*`alpha25_1890'
						- _Iyear_1900*`alpha25_1900'
						- `theta25'*ln_tpop
				if (banks==0);
replace eta0ct25 = ln(`bestdivide25'/(2*tpop))-`alpha25' 
						- _Iyear_1880*`alpha25_1880'
						- _Iyear_1890*`alpha25_1890'
						- _Iyear_1900*`alpha25_1900'
						- `theta25'*ln_tpop
				if (banks1902==0) & year == 1900;				
by gisjoin1890: egen etai_c25 = mean(eta0ct25);
drop eta0ct; 

gen bestdivide_pc = `bestdivide'/tpop;
gen bestdivide25_pc = `bestdivide25'/tpop;
gen bestdivide50_pc = `bestdivide50'/tpop;

sort gisjoin1890 year;
tempfile tempdata;
save `tempdata', replace;

/*****Make Multiple imputations ***/
tempfile temp2;
foreach do25 in "" "25"  {;/* repeat for using 1900 data, and for replacing 1900 with 1902*/

use `tempdata', clear;
if "`do25'" == "25" {;
	foreach var in lowcapital highcapital needscapital etai_c {;
		replace `var' = `var'25;
	};
};

/*Reduce data to only rows, columns necessary*/
keep gisjoin1890 year lowcapital highcapital /* banks50c banks25c banks*/ needscapital etai_c ;
keep if needscapital;
reshape wide /*banks banks50c banks25c*/ needscapital lowcapital highcapital, i(gisjoin1890) j(year);
local stubnames;
forvalues ii = 1/`replications' {;
		quietly gen eta`ii'_c = .;
		local stubnames `stubnames' eps`ii'_c;
};

foreach thisyear in 1870 1880 1890 1900 {;
	forvalues ii = 1/`replications' {;
		quietly gen eps`ii'_c`thisyear' = .;
	};

};
/*Call mata program, which accesses stata data*/	
mata makedraws(`Ninitial');


drop /* banks* banks50c*/ needscapital* lowcapital* highcapital* etai_c;
reshape long `stubnames', i(gisjoin1890) j(year);
sort gisjoin year;

save `temp2', replace;
use `tempdata', clear;
merge 1:1 gisjoin1890 year using `temp2';
drop _merge;


/*Create base value for excess capital that is 0 for those that do not have excess capital, and missing for those that do. This variable will be imputed*/

gen excess_capital_pc =0 if needscapital==0;
gen opt_capital_pc = capitalstock_pc if needscapital==0;


mi set wide;
mi set M=`replications';
mi register regular  gisjoin1890- eps3_c;
mi register imputed excess_capital_pc opt_capital_pc;

forvalues ii = 1/`replications' {;
	replace _`ii'_excess_capital_pc = excess_capital_pc;
	replace _`ii'_excess_capital_pc = capitalstock_pc 
					- exp(`alpha`do25''+ eta`ii'_c + eps`ii'_c 
						 	+ _Iyear_1880*`alpha`do25'_1880'
							+ _Iyear_1890*`alpha`do25'_1890'
							+ _Iyear_1900*`alpha`do25'_1900'
							+ `theta`do25''*ln_tpop ) if needscapital;
	
		
	replace _`ii'_opt_capital_pc = opt_capital_pc;
	replace _`ii'_opt_capital_pc =  exp(`alpha`do25''+ eta`ii'_c + eps`ii'_c 
						 	+ _Iyear_1880*`alpha`do25'_1880'
							+ _Iyear_1890*`alpha`do25'_1890'
							+ _Iyear_1900*`alpha`do25'_1900'
							+ `theta`do25''*ln_tpop ) if needscapital;					 
	
};
/*Save MI file*/
compress;
cd `OUTDIR';
save MI_National_banks`do25'_`sampleselect'_v0, replace;
/*** Naming convention 
v0 allows scale effects
v3 does not allow scale effects in population (theta set to 0, no ln_tpop in MLE)
***/		
};/* End loop for creating both banks with 25 capital and 50 capital in 1900*/
/********** End create MI file  ****/

log close;
exit;

/*****Draw from normal repeatedly until get within bounds******/
/*Old version, some of the bounds are so tight that getting a draw in bounds is very unlikely						
gen eps_ct =.;
gen eta_c =.;
local anotherdraw =1;
local xx = 1;
/*pause on;*/
quietly {;
while (`anotherdraw' ==1) {;
	/*Draw epsilons for all counties and years that still need a valid draw*/
	replace eps_ct = rnormal(0,`sigma_eps') if stillneedscapital;
	/*Draw eta's for all counties that still need a valid draw eta's are the same for all years in a county*/
	by gisjoin1890: replace eta_c = rnormal(0,`sigma_eta') if _n==1 & stillneedscapital;
	by gisjoin1890: replace eta_c = eta_c[1] if stillneedscapital;
	
	replace needscapital = (insample & (year >=1870 & year <=1900) &
					(/*Counties with a 50 bank*/ (banks50c>0 & banks50c<.) & 
						( (lowcapital > eps_ct + eta_c) | (highcapital<eps_ct + eta_c)) 
					)	| (/*Counties without a bank*/ (banks==0) & 
						 (highcapital < eps_ct + eta_c) )
					); 
	drop stillneedscapital;
	by gisjoin1890: egen stillneedscapital = max(needscapital);
	sum stillneedscapital;
	noisily display "Iteration `xx': `r(sum)'";
	local anotherdraw = r(max);
	/*if `xx' >= 3000 {;
		pause;
	};*/
	local xx = `xx'+1;
	
}; };
*/
